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Abstract 


We have developed a novel method for co-adding multiple under-sampled images that combines the iteratively 
reweighted least squares and divide-and-conquer algorithms. Our approach not only allows for the anti-aliasing of 
the images but also enables Point-Spread Function (PSF) deconvolution, resulting in enhanced restoration of 
extended sources, the highest peak signal-to-noise ratio, and reduced ringing artefacts. To test our method, we 
conducted numerical simulations that replicated observation runs of the China Space Station Telescope/ the VLT 
Survey Telescope (VST) and compared our results to those obtained using previous algorithms. The simulation 
showed that our method outperforms previous approaches in several ways, such as restoring the profile of extended 
sources and minimizing ringing artefacts. Additionally, because our method relies on the inherent advantages of 
least squares fitting, it is more versatile and does not depend on the local uniformity hypothesis for the PSF. 
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However, the new method consumes much more computation than the other approaches. 
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1. Introduction 


As technology has advanced, computer software has 
become increasingly important in the field of image proces- 
sing. To simplify data processing, an Image Acquisition 
System (IAS) is often utilized for image digitization. 
Essentially, an IAS acts as a digitizer that converts continuous 
signals into digital ones by recording data with detectors, such 
as pixels, which appear in a mosaic-like pattern known as 
pixelation. However, due to the diffraction limit of the optics 
equipment, the images captured by the IAS are band-limited. 
This means that the signal recorded by the image has a 
maximum spatial frequency or resolution, as described in 
Fruchter (2011). A band-limited signal can be fully recon- 
structed by high-sampling imaging, according to the Nyquist 
sampling theorem. Nevertheless, many astronomical images 
do not meet the Nyquist sampling criterion due to technical or 
economic limitations, resulting in under-sampled images 
(Fruchter & Hook 2002). 

Under-sampled images suffer from an aliasing effect, which 
blurs all signals below a sampling interval. To achieve a 
resolving power that approaches the diffraction limit, the 
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sampling rate needs to be increased. This can be accomplished 
by combining multiple exposures using various coaddition 
methods, such as shift-and-add (Bates & Cady 1980; Farsiu 
et al. 2004b), Drizzle (Fruchter & Hook 2002), Super-Drizzle 
(Takeda et al. 2006), IMCOM (Rowe et al. 2011), iDrizzle 
(Fruchter 2011), SPRITE (Ngolé Mboula et al. 2015), and 
fiDrizzle (Wang & Li 2017), as well as iterative back-projection 
(IBP, Irani & Peleg 1993; Symons et al. 2021). Drizzle has 
become the standard for combining images taken by the 
Hubble Space Telescope (HST) and the James Webb Space 
Telescope (SWST). Some Drizzle-based methods are widely 
used to restore fine details of under-sampled multi-exposures 
and fuse images from different equipment. SPRITE and IBP are 
developed to reconstruct the Point-Spread Function (PSF) 
using stars in even one exposure. 

Due to the diffraction of light, the image of a point-like source 
in the image plane is not actually a point. Instead, it appears as an 
extended object, known as the PSF effect. This degradation of the 
observed image is caused by the convolved PSF, but PSF 
deconvolution technology can improve resolution by compensat- 
ing for this numerically. Along with resolution improvement, 
PSF deconvolution also enhances contrast and reduces noise 
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(Sage et al. 2017). There are various algorithms for PSF 
deconvolution (Starck et al. 2002), including the Fourier-quotient 
method, CLEAN method (Hégbom & Cornwell 1974), Bayesian 
approach [including Landweber method (Landweber 1951), 
Richardson—Lucy algorithm (Richardson 1972; Lucy 1974; 
Shepp & Vardi 1982)], wavelet-based deconvolution, and 
super-resolution techniques (Gerchberg 1974; Hunt 1994; Elad 
& Feuer 1999; Lauer 1999; Capel & Zisserman 2003; Park et al. 
2003; Farsiu et al. 2004a; van Ouwerkerk 2006; Tian & 
Ma 2011; Nasrollahi & Moeslund 2014; Yue et al. 2016; Symons 
et al. 2021). The PSF deconvolution has numerous benefits, such 
as its applicability to even the simplest optical setup, reduction of 
financial costs, and streamlining of the acquisition pipeline. 
However, when noise contaminates the PSF-convolved image, 
PSF deconvolution becomes an ill-posed problem (Starck et al. 
2002; Sage et al. 2017). Regularized methods are often employed 
to generate an approximation (Ng et al. 2007; Takeda et al. 
2007, 2009; Yuan et al. 2010; Babacan et al. 2011; Su et al. 
2012; Zhang et al. 2012; Liu & Sun 2014), which is effective and 
flexible in reducing noise amplification, ringing effect, and flux 
divergence. 

However, most previous works aim to either the PSF 
deconvolution for single exposure (e.g., CLEAN, Landweber, 
Richardson—Lucy, maximum entropy algorithm, etc.) or to 
multiple exposures coaddition but without PSF deconvolution 
(e.g., shift-and-add and Drizzle based methods). Only a few 
works mention the under-sampled multi-exposures coaddition 
(MEC) with PSF deconvolution, e.g., Stark-Pantin: Equation 
(53) in Starck et al. (2002) and UPDC: Equation (10) in Wang 
et al. (2022). In this paper, we systematically study how to 
achieve superresolution in the iterative MEC by anti-aliasing 
and PSF deconvolution coaddition (AAPD). We test several 
recovery methods with numerical simulations. 


2. Imaging Model and the Least Squares 


The observation image records not only the shape or flux from 
the objects of study but also a set of combined observational 
effects, e.g., blurring effects (PSF, seeing, vignetting, etc.), 
sampling effect (pixelation, image field distortion), noise effect 
(equipment noise, environment noise), etc. Therefore, a reason- 
able imaging model should take those observational effects into 
account. 

Let B be an image degrading operator (i.e., blurring) which 
represents a combination of a series of image operations e.g., 
PSF convolution, pixelation, etc. The formation of an under- 
sampled observation image G with the size of U x V pixels can 
be formally expressed as 


G = N{B{O}}, (1) 


where © is the intrinsic, continuous surface brightness of the 
sky. Assuming that © contains a fine grid of size M x N, we can 
represent the downgraded effects of the image, except for the 
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Table 1 

Important Variables Used in the Text 
Variable Description 
oO intrinsic, continuous surface brightness of the sky 
Gk the kth exposure 
T fully sampled high resolution (or definition) image 
Pk PSF measured from the kth frame 
W weight for flux allocation 
G sample for pixels on the same WCS position 
B down-sampling factor 
0 noise at the pixel 
N{-} noise operator 
B{-} blurring operator 


Pgri} component-wise projection operator for constraints set 


noise, using B {OQ}. The noise contamination on the downgraded 
image is described by N. To recover the original image O from 
G, we need to know the combined effects of N{B{-}} in 
Equation (1) in advance. In this study, we focus on the recovery 
method and assume that the PSFs and the position shift on all 
dithered exposures are well measured beforehand. Table 1 
provides a list of important symbols and their representations. 

Assuming that we have L blurred (PSF convolved and/or 
under-sampled) exposures {Gj, Go... GL}? with Ux V pixels 
for the same object and the known blurring effects B{-}, can we 
mock another set of frames {Fi, F2... FL} to approach the 
observations, then constrain the original image? To answer this 
question, we try to use the least squares to calculate a x? for all 
pixels between the observation and mock samples. Theoreti- 
cally, the x? can be represented as 


3 N (Fia Gan 


2 
has 


: (2) 


k,ij 


where F = B{T} and T is a target image with size of M x N 
fine grids, i.e., the estimation of the original image O, while 
Orij is the noise at the pixel of the ith row and jth column of the 
kth observational frame. Note that the summation of the 
squares is over L x U x V pixels. Since the blurring effects 
{-} are known, minimizing the xX? results in an estimation of 
the target image 7 . By differentiating x? with respect to 7 and 
setting the derivative to zero, we have 
Oy? _ 
OT 
It is very difficult to solve Equation (3) directly, especially 
for huge amounts of pixels. Actually, Equation (3) is a system 
of linear equations with Mx N unknowns and Lx Ux V 
conditions. As a kind of fallback solution, we employ a divide- 


and-conquer algorithm to solve Equation (3) on each target grid 
(m, n) at a time. The Drizzled grids are used as the initial 


(3) 


12 Frame Gy (k = 1...L) is obtained from the kth observation. 
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Figure 1. Resampling mechanism. The target grid TO, is shown as the fine (in black), the observation pixels G,;,; are shown as coarse (in red). The flux weight 
W.i.j,»,q (in pink) is determined by the clipped area between the observation pixels and the target grid. Here the TS space (p, q) includes all target cells that overlap 


with the observation pixel Gz,;j (blue region). 


values. The target grid (m, n) will be renewed after solving 
Equation (3). The program checks if the output meets the preset 
threshold'* when all grids are updated. If not, the updated grids 
will take the place of the initial values in Equation (3). Then it 
repeats the above steps until the updated grids meet the 
threshold. Therefore, Equation (3) is solved in an iterative 
approach as follows: 


3 02 g È Maio Penna 


kij Onis pa St 
TS 


x 
a) Pk.p,q,m,n) D Wer, ijp,aPkp.q,m,n 
Pq 


pirn- 


Os Gr kij 
= tea eae: (4) 
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where P; is a fine PSF measured from the kth frame. The PSF 
must have the same sampling rate as the target image. Wki j,p,q 
is a flux allocation weight that stands for how many areas of an 


13 Th this work, we suggest using the x° as the threshold. 


observation pixel (k, i, j) are overlapped by the target pixel 
(p, q), shown as Figure 1. In this work, the overlapped polygon 
of two grids is clipped using a polygon clipping algorithm 
Clipper2.'* The area of the clipped polygon is calculated using 
the Green formula. Indexes r and r+ 1 represent the sequence 
number of the iteration step. OS is a space that includes all 
observation pixels that are overlapped with PSF P, which is 
centered at the target pixel (m, n), while TS includes all target 
pixels that are a with the OS observation pixels. 
Let Y, ij,mn =r q Wk, LP, qFx, P.qmn and oP = 


k,i,j,m,n 
D q Wiij.p, PD r Pepis TO ;, Equation (4) can be rewritten as 


OS %kij,m,n (r) 
eae kij gp, Dris — Pkijmn) 
tä fig 
Tin = Tni T v, (5) 
jimin 
D sid OR, 


which is an iterative solution to the least squares. po) jum 18 the 


mimic flux of pixels that are overlapped by the target grid T“”? 


mn* 
In order to eliminate contaminations from cosmic rays or 


14 https: //github.com/AngusJohnson /Clipper2 
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Figure 2. Cosmic-ray removal. For one of the exposures, there are four circles in the top panels and bottom left. The bottom right panel shows all the mismatched 
cases (including 50 exposures). Circles in green stand for the misidentification cases (They are not real cosmic rays), while blue for the omissions. 


abnormal pixels, the weight Ox ;j is dynamically adjusted with 
the iterations, which results in robust regression of the 
algorithm. This is actually an application of the Iteratively 
Reweighted Least Squares (IRLS) method in image stacking. 
We call this divide-and-conquer IRLS method i.e., Equation (5) 
as Dirles. The initial value of the iteration 7® is set as the 
Drizzled real exposures to automatically introduce telescope 
effects such as vignetting, instrument noise, and saturation 
overflow. Therefore, these effects do not need to be simulated 
anymore in the Dirles. 

In literature, there are two algorithms besides the Dirles that 
can achieve the anti-aliasing and PSF deconvolution coaddi- 
tion simultaneously. These algorithms are called Starck—Pantin 
(Starck et al. 2002) and UPDC (Wang et al. 2022). Both 
Starck—Pantin and UPDC are forward modeling methods that 
infer the original image by maximizing the likelihood between 
observations and mocks. The algorithms like Richardson— 
Lucy, Starck—Pantin, and UPDC use FFT to convolve the PSF. 
However, the Dirles algorithm convolves the PSF in the real 
space, making it more suitable for situations where the PSF 
changes dramatically. Although it consumes much more 
computation than the previous algorithms, it is a better choice 
for such situations. 

Drizzling the multi-exposures of an ideal point source can 
form a blurring effect in the target image (grid), which 
originates from the pixelation. It is called pixelation blur to 
distinguish it from the ordinary PSF (see Wang et al. 2022 for 
more details). The pixelation blur is not homogeneous, even 


not continuously changing in the whole target image. Because 
different point sources have different positions therefore 
different cases of pixels coaddition. Due to the heterogeneity 
of pixelation blur, theoretically, it is not a good choice to 
deconvolve the PSF by using ordinary methods e.g., Land- 
weber or Richardson—Lucy algorithm on the Drizzled image. 


3. Image Completion: Cosmic Rays and Bad Pixels 
Replacement 


Instead of using the Drizzle method for cosmic-ray removal 
(Fruchter & Hook 2002), we apply an improved statistical 
algorithm that identifies and replaces abnormal regions with 
reasonable values. The new algorithm can significantly reduce 
the mismatched cosmic rays (~50%) in our China Space 
Station Telescope - Multi-channel Imager (CSST-MCT) mock 
test, e.g., Figure 2. The process of image completion involves 
the following steps. 


1. To start, let us create a target grid, denoted by T,,,,, with a 
1:1 sampling rate. This grid will be utilized to resample 
all L exposures Gi; (kK=1...L) based on the WCS 
parameters specified in their headers. 

2. For each target pixel located at position (m, n), we collect 
a statistical sample of overlapping exposure pixels, 
denoted by Go(m,n). We gather a total of Mx N 
samples. These statistical samples are then utilized to 
calculate the median value, providing the first estimation 
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of an image that is free of cosmic rays or bad pixels, 
known as Tmia(m, n). 

3. The mean, uo(m, n), and standard deviation, oo(m, n), are 
also calculated for each sample Go(m, n). We use these 
values to identify cosmic-ray candidates by checking if 
the exposure pixels with a flux of Gri; > Tmia(m, n) + 
Vo X (m, n) meet the preliminary voting condition. 
Here, v is an artificial threshold set to 5 in our work. If a 
region occasionally has multiple cosmic rays, the mean 
Holm, n) can be overestimated. Therefore, we replace the 
mean pig(m, n) with the median Zmia(m, n) in the 
preliminary voting condition. 

4. Typically, an exposure pixel, Gk, j, is covered by 1 to 5 
neighboring target pixels that act as voters” in 
determining whether a cosmic-ray candidate is present. 
If a cosmic-ray candidate receives all the votes of the 
voters in the district, it is flagged as a cosmic ray. We 
remove the cosmic rays from the sample Go(m, n) to 
obtain a relatively pure sample, denoted by G,(m, n). 

5. Once again, we calculate the mean, uı(m, n), and 
standard deviation, o,(m, n), for each sample G,(m, n). 
We then perform a series of similar operations as 
described above, but with a refreshed condition 
Gk ij > (m, n) + ni x alm, n), where v is set to 5. 
This allows us to remove cosmic rays from sample 
Gı(m, n), resulting in a reduced sample called G2(m, n). 
Please note that some exposure pixels with flux 
Grij > (m, n) + vo X a(m, n) are removed for their 
cosmic-ray identification, but those with flux Grij< 
u(m, n) — vo X oj(m, n) are not. This means that using 
the mean (m, n) and standard deviation o2(m, n) taken 
from sample G2(m, n) to generate a Gaussian distribution 
is biased. To avoid that issue, we fit the sample G2(m, n) 
to a Gaussian profile using the nonlinear least squares 
fitting code MPFIT."° 

6. Each bad pixel or cosmic ray exposure pixel will be 
replaced with a random value generated by its corresp- 
onding Gaussian profile. 


The above method is not valid for pixels that are overlapped by 
only one or two exposures due to the lack of a statistical 
sample. These cases usually occur at the marginal area of the 
target image Tmn, and can be ignored. Figure 2 shows a case 
for cosmic-ray removal. The circles and ellipticals denote the 
mismatched cosmic rays. The bottom right panel shows the 
total number of mismatched cases for 50 exposures. The blue 
circles and ellipticals represent cosmic rays that were omitted 
due to their low flux blending with the background, while the 
green circles indicate cases of mistaken removal caused by the 
threshold set-up, i.e., a small probability event below the 
threshold. 


15 https: //pages.physics.wisc.edu /~craigm/idl/cmpfit.html 
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4. Simulation Tests Based on the CSST-MCI Mock 
Pipeline 

By employing the CSST-MCI mock pipeline, we take one 
HST Wide Field Camera 3 (WFC3) stamp from a real galaxy 
database of Galsim (Rowe et al. 2015) as our mock input (or 
original image, ground truth), to generate mock exposures with 
random shifts and rotations for the g-band. In the CSST-MCT 
mock, four blurring effects are taken into account: PSF 
convolution, cosmic ray, down-sampling, and Gaussian/ 
Poisson noise. Following Wang et al. (2022), a non-strict 
positivity constraint is adopted: 


THY) THD SO 


7 (6) 
TÒ THD<O 


PeT) = { 


where PgH{TC+D} is a component-wise projection of T+) 
onto the set R* (not strictly). For the convenience of 
comparison, we assumed that a PSF is uniform in a local 
region. The region size is determined by the spatial growth rate 
of the PSF. It is a local uniformity hypothesis for the PSFs. 


4.1. Mock-I: Faint Sources Detection 


This case centers on the application of mock and multi- 
exposure coaddition to CSST-MCT, i.e., Figure 3. The image in 
the lower left panel serves as the ground-truth input, extracted 
from an HST observation and showcasing no less than 15 
substructures'® on the central galaxy. The upper left panel features 
one of the g-band exposures (512 x 512 pixels, with a down- 
sampling factor G=2, convolved PSF size 128 x 128 pixels), 
which displays two easily identifiable substructures. The Drizzled 
image for 50 exposures is displayed in the upper right panel, 
revealing roughly seven substructures. Lastly, the final panel 
showcases the coadded image generated by this work, where all 
15 substructures on the central galaxy are clearly visible. There 
are seven stars in the field of view, corresponding to seven 
pixels with high flux in the ground-truth input. The stars are 
visible in the rest panels but with different full half-maximum 
widths (FHWM). 


4.2. Mock-II: Strong Lensing Mock and Images 
Coaddition 


Figure 4 shows a strong gravitational lensing mock and 
recoveries from four kinds of coaddition methods. The panel in 
the bottom left displays the ground-truth image generated by a 
lensing model and used as the mock input. Notably, the strong 
lensed giant arc within the dashed annulus is known as the 
Einstein ring. Moving to the top middle panel, we see an 
example of one of the CSST-MCI mocked exposures in the 
r-band with a down-sampling factor G=4. Finally, the 


16 The substructures of galaxies are detected by the software SourceXtractor++ 
(https: / /sourcextractorplusplus.readthedocs.io/en /latest/Introduction.html). It may 
be point-like sources or spiral arms in the figure. 
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15 substructures 
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Figure 3. CSST-MCI mock and multi-exposure coaddtion. The bottom left panel shows the ground-truth image as the mock input, which is extracted from the HST 
observation with at least 15 substructures on the central galaxy. One of the g-band exposures with two obvious substructures is illustrated on the top left. There are 
about seven substructures on the top right panel (the Drizzled image for 50 exposures). All the 15 substructures on the central galaxy can be found on the last panel, 


which is coadded by this work. 


remaining panels showcase four distinct coaddition methods 
that were used to recover images from the exposures 
(128 x 128 pixels, convolved PSF size 128 x 128 pixels). 
These methods are referred to as Drizzle (top right), UPDC 
(bottom left), Starck-Pantin (bottom middle), and Dirles 
(bottom right). Visually, UPDC and Dirles manifest better 
results than the Drizzle or the Starck—Pantin in the Einstein 
ring recovery from 100 mocked exposures. 

Quantitatively, we measure the peak signal-to-noise ratio 
(PSNR) within the region of the dashed annulus in Figure 4 for 
each iteration of different coaddition methods. The result is 
shown in Figure 5. The methods tested are Drizzle (with a 
PSNR of 22.16), UPDC (in Red), Richardson—Lucy (in Green), 
Starck—Pantin (in Blue) and Dirles (in Black). The curves end 
at the optimal iteration with the highest PSNR. The Dirles 
achieves the highest PSNR among all the recoveries with less 
than 40 iterations. Richardson—Lucy method consumes the 
most iterations but the least computation (see Table 2). Note 


that the same PSFs are used in the PSF deconvolutions as those 
for the mock input. 


4.3. Mock-III: VST Mock and Image Restoration 


To test the recovery methods, we have mocked the 100 
frames r-band VOICE — CDFS — 1 multi-exposures data from 
the OmegaCAM of the Very Large Telescope Survey 
Telescope in the European Southern Observatory. To reduce 
computation time, we extracted a square region that is centered 
at [R.A. = 53°1511, decl. = —27°7175] and has a size of 
101 x 101 observation pixels. There is a high-resolution image 
(goodss_3dhst_F606W_sci.fits) from the HST in this region. 
The high-resolution HST image has been binned 
1.77858 x 1.778 58 times along the two-dimensions. 

In our simulation process, we use a group of Moffat profiles 
to replicate extended sources such as galaxies. To mimic the 
central point source (the star) on the target grid, we adopt a 
single pixel. By adjusting the Moffat profiles’ parameters and 
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CSST-MCI reband mock 
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Figure 4. A strong gravitational lensing mock. The bottom left panel shows the ground-truth image as the mock input, which is generated by a lensing model. One of 
the CSST-MCI mocked exposures in the r-band is illustrated on the top middle panel (with a down-sampling factor 8 = 4). Four kinds of coaddition methods 
recovered images from exposures are shown on the rest panels: Drizzle (top right), UPDC (bottom left), Starck—Pantin (bottom middle) and Dirles (bottom right). 


the central point’s brightness, and convolving a PSF from the 
binned HST image on the same sky field, we can fit a target 
image to the binned HST image. The target image is shown on 
the top left panel of Figure 6, and we consider it to be true. 

To produce the mock samples, we first combine the fine 
PSFs and true PSFs measured in the 100 VOICE frames. The 
fine PSFs, which have a dimension of 59 x 59 pixels, are 
generated by a Principal Component Analysis (PCA) based 
method developed by our team’s cooperators Nie et al. 
(2021a, 2021b). Next, we shift the image using the same 
position shifts as the VOICE frames and down-sample the true 
image to the same coarse grids as in the VOICE frames. 
Finally, we add the Poisson and Gaussian noise from each 
VOICE frame to the corresponding mock sample, as shown in 
the middle panel of the top array of Figure 6. 

To make the results more realistic, we create sets of 
realizations by modifying the random seed in the Poisson and 
Gaussian processes. For each set of realizations, we reconstruct 
the mock samples using the methods mentioned in the text, 
namely Drizzle, UPDC, Starck—Pantin, and Dirles, up to the 
optimal number of iterations!” in the 2 x 2 finer target grid. 

According to the results depicted in Figure 6, the Dirles 
approach appears to yield superior reconstructions in compar- 
ison to earlier methods. This can be attributed to its ability to 
effectively recover object shapes, minimize the ringing effect, 


17 The iteration has the highest PSNR. 


and suppress background noise. It should be noted, however, 
that Dirles may disregard certain low SNR regions, such as the 
tails of galaxies. To circumvent this potential concern, 
decreasing the number of iterations and gradually balancing 
noise reduction with signal enhancement is advisable. 

Weak gravitational lensing is primarily concerned with 
measuring the shear. This involves measuring the deformation 
of the background galaxies compared to the randomly aligned 
ones. By doing so, it is possible to constrain the properties of 
the foreground lensing objects. Therefore, the ellipticity of 
galaxies is the signal that researchers are interested in 
extracting. Following Hirata & Seljak (2003), the ellipticity 
of an object is defined as 


e4 = (Mx. — Myy)/(Mex F Myy) 
ex = 2M; ,/ (Mx F Myy) (7) 


where M,; represents the moments, the spin-2 tensor e = (e4, 
ex) is the so-called ellipticity tensor. To address the problem of 
divergence, we apply a circular Gaussian weighting function to 
the largest galaxy, with a weight radius of 7,,. The variance of 
shape parameters is then plotted in Figure 7 against the weight 
radius r„. By analyzing the visual illustration in Figure 6 and 
the quantitative results in Figure 7, we can conclude that the 
Dirles and UPDC methods outperform other methods in 
reconstructing shape parameters for extended sources. It is 
worth noting that the Drizzle method produced the largest 
deviation from the ground truth. 
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Figure 5. A comparison of the PSNR from different coaddition methods: Drizzle (PSNR = 22.16), UPDC (in Red), Richardson—Lucy (in Green), Starck—Pantin (in 
Blue) and Dirles (in Black). The curves end at the optimal iteration with the highest PSNR. 


Table 2 
Computational Resources Consumed by Various Coaddition Approaches During Each Iteration 
Simu No. Drizzle Richardson—Lucy Starck—Pantin UPDC Dirles 
Mock-I (512 x 512 x 50, 6 = 2) 3.7 s 220.4 ms 7.9 s 7.8 s 272.2 s 
Mock-II (128 x 128 x 100, 8 = 4) 1.9 s 90.1 ms 3.8 s 3.9 s 144.3 s 
Mock-MI (101 x 101 x 100, 8 = 2) 0.3 s 0.6 s 0.6s 5.458 


In the field of photometry, the reconstructed source profile 
plays a pivotal role in estimating the performance of recovery 
algorithms. To analyze the profiles for both the central star and 
the largest nearby galaxy, we have presented graphical 
representations in Figures 8 and 9, respectively. Our analysis 
based on Figures 6, 8, and 9 demonstrates that Dirles 
outperforms other methods by effectively recovering high peak 
intensities in point sources, providing superior profiles for 
extended sources, and minimizing undesirable ringing arte- 
facts. It is noteworthy that, for a given recovery method, the 
optimal number of iterations exhibits minimal variation across 
different realizations. 

We have conducted a comprehensive study of the computa- 
tion time required by different coaddition methods during each 
iteration. The results are presented in Table 2. All calculations 
were performed on a two-socket AMD server equipped with 


two CPUs of EPYC 7763 @2.45-3.5 GHz and 1 TB DDR4 
memory. 

The Richardson—Lucy method has the lowest computation 
cost as it mainly executes the FFT in the calculation, which 
increases with an O(N logN) trend. On the other hand, 
Drizzle consumes the most time on the exposure up-sampling 
process. Both Starck—Pantin and UPDC involve similar 
operations such as PSF convolving (via FFT), down- 
sampling, comparing, and up-sampling. Hence, they have 
almost the same computation cost, which is mainly con- 
tributed by the sampling. 

Dirles has the highest computation cost compared to the 
other methods. Instead of the FFT-based convolution, its PSF 
convolution is performed in real space, which significantly 
increases computational time. Apart from the size of input 
exposures, the computation also depends on the dimension of 
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Figure 6. Recoveries from the VST mock. The top left panel displays the ground truth, which is a model fitting to the binned HST observation. One of the VST mock 


exposures is shown on the top middle panel. The top right panel displays the image recovered by using Drizzle method. The three panels at the bottom, from left to 
right, represent the recoveries by UPDC, Starck—Pantin, and Dirles methods respectively. 
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Figure 7. Shape parameters comparison between the true and the recoveries: left for e, and right for ex. 


the PSF provided in advance. However, the Dirles method 
stands out from other approaches as it does not assume local 
uniformity for the PSF. Consequently, it is a better fit for 
scenarios where the PSF undergoes significant changes within 


an exposure space. With Dirles, we can convolve a PSF map 
that is different everywhere in the field of view. D. Liu et al. 
(2023, in preparation) are constructing such a PSF map by 
using a set of basic functions for fitting. 
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Figure 8. The left panel shows the normalized profiles of the central star drawn in different colors to differentiate between them: Drizzle (in Magenta), UPDC (in Red), 
Richardson—Lucy (in Green), Starck—Pantin (in Blue), and Dirles (in Black). The right panel is a zoomed-in version that focuses on the area around the zero flux. This 


view highlights the ringing effect and other details around the point source. 
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Figure 9. Profiles of the largest galaxy from ground truth (in Cyan) and five recoveries: Drizzle (in Magenta), UPDC (in Red), Richardson—Lucy (in Green), Starck— 


Pantin (in Blue) and Dirles (in Black). 


5. Conclusions 


In this article, we propose a novel AAPD method, named 
Dirles, which is based on the least squares fitting technique. 
Following a series of rigorous simulation tests conducted on 
these algorithms, our findings can be summarized as follows. 


1. In detecting faint sources, Dirles recovered all 15 
substructures while Drizzle missed 8. 
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2. The Dirles approach outperforms previous methods by 
achieving the highest PSNR, recovering object shapes, 
minimizing the ringing effect and suppressing noise. 
However, it may ignore low SNR regions. To address 
this, gradually reduce iterations and balance noise 
reduction with signal enhancement. 

3. The Dirles method outperforms other methods in 
reconstructing shape parameters for extended sources. 
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In radiometry, the Dirles method performs the best in 
terms of recovering the highest peak in point sources, 
providing the best profile for extended sources, and 
reducing the ringing effect. 

4. Compared to the other methods, Dirles has the highest 
computation cost. However, it stands out from other 
approaches as it does not assume local uniformity for the 
PSF, making it a better fit for scenarios where the PSF 
undergoes significant changes within an exposure space. 


After extensive simulations, it can be concluded that Dirles 
outperforms previous works in restoring point/extended 
sources, maintaining extended source shapes, and reducing 
ringing despite high computation consumption. 

In the future, numerous upcoming telescopes will commence 
astronomical observations, such as NASA’s Wide Field 
Infrared Survey Telescope (WFIRST), the European Space 
Agency’s Euclid mission, the National Science Foundation- 
funded Large Synoptic Survey Telescope (LSST), and China’s 
Space Station Optical Telescope (CSST). These advanced 
instruments will generate a vast volume of imaging data. 
Effectively processing these images while maintaining high 
fidelity is an imminent necessity. As the implementation of 
Dirles algorithm uses a divide-and-conquer approach, it may be 
possible to accelerate it to a considerable extent using GPU- 
based parallel computing. We are confident that there is a 
significant potential for improving the speed and effectiveness 
of the Dirles method even further. 
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In this study, a cluster is used with the SIMT accelerator 
made in China. The cluster includes many nodes each 
containing two CPUs and four accelerators. The accelerator 
adopts a GPU-like architecture consisting of a 16GB HBM2 
device memory and many compute units. Accelerators 
connected to CPUs with PCI-E, the peak bandwidth of the 
data transcription between main memory and device memory is 
16 GB s™'. 
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